function [lambda,weight] = quadpts3(order)
%% QUADPTS3 quadrature points in 3-D.
%
% [lambda,weight] = quadpts3(order) return quadrature points in a
% tetrahedron with given order (up to 5) in the barycentric coordinates.
%
% The output lambda is a matrix of size nQ by 3, where nQ is the number of
% quadrature points. lambda(i,:) is three barycentric coordinate of the
% i-th quadrature point and lambda(:,j) is the j-th barycentric coordinate
% of all quadrature points. The x-y-z coordinate of the p-th quadrature point
% can be computed as 
%
%    pxyz = lambda(p,1)*node(elem(:,1),:) ...
%         + lambda(p,2)*node(elem(:,2),:) ... 
%         + lambda(p,3)*node(elem(:,3),:) ... 
%         + lambda(p,4)*node(elem(:,4),:);
%
% The weight of p-th quadrature point is given by weight(p). See
% verifyquadpts for the usage of qudrature rules to compute integrals over
% triangles.
% 
% References: 
%
% * Jinyun, Y. Symmetric Gaussian quadrature formulae for tetrahedronal
% regions. Comput. Methods Appl. Mech. Engrg.. 43(3):349--353, 1984.
%  
% See also quadpts1, quadpts, verifyquadpts
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details. 

if order>5
    order = 5;
end
switch order
    case 1    % Order 1, nQuad 1
        lambda = [1/4, 1/4, 1/4, 1/4];
        weight = 1;
    case 2    % Order 2, nQuad 4
        alpha = 0.5854101966249685; 
        beta =  0.138196601125015;
        lambda = [alpha beta beta beta; ....
                  beta alpha beta beta; ...
                  beta beta alpha beta; ...
                  beta beta beta alpha];
        weight = [1/4, 1/4, 1/4, 1/4];
    case 3    % Order 3, nQuad 5
        lambda = [1/4 1/4 1/4 1/4; ...
                  1/2 1/6 1/6 1/6; ...
                  1/6 1/2 1/6 1/6; ...
                  1/6 1/6 1/2 1/6; ...
                  1/6 1/6 1/6 1/2];
        weight = [-4/5, 9/20, 9/20, 9/20, 9/20];
    case 4    % Order 4, nQuad 16
        alpha1 = 0.7716429020672371; 
        beta1 =  0.7611903264425430e-1;
        w1 = 0.5037379410012282e-1;
        alpha = 0.4042339134672644;
        beta = 0.7183164526766925e-1;
        gamma = 0.11970052777978019;
        w2 = 0.6654206863329239e-1;
        lambda = [alpha1 beta1 beta1 beta1; ....
                  beta1 alpha1 beta1 beta1; ...
                  beta1 beta1 alpha1 beta1; ...
                  beta1 beta1 beta1 alpha1; ...
                  alpha alpha beta gamma; ...
                  alpha alpha gamma beta; ...
                  alpha beta alpha gamma; ...
                  alpha beta gamma alpha; ...
                  alpha gamma beta alpha; ...
                  alpha gamma alpha beta; ...
                  beta alpha alpha gamma; ...
                  beta alpha gamma alpha; ...
                  beta gamma alpha alpha; ...
                  gamma alpha alpha beta; ...
                  gamma alpha beta alpha; ...
                  gamma beta alpha alpha];                  
        weight = [w1, w1, w1, w1, ...
                  w2, w2, w2, w2, w2, w2, ...
                  w2, w2, w2, w2, w2, w2];
    case 5    % Order 5, nQuad 17
        alpha1 = 0.7316369079576180; 
        beta1 =  0.8945436401412733e-1;
        w1 = 0.6703858372604275e-1;
        alpha = 0.4214394310662522;
        beta = 0.2454003792903000e-1;
        gamma = 0.1325810999384657;
        w2 = 0.4528559236327399e-1;
        lambda = [1/4, 1/4, 1/4, 1/4; ...
                  alpha1 beta1 beta1 beta1; ....
                  beta1 alpha1 beta1 beta1; ...
                  beta1 beta1 alpha1 beta1; ...
                  beta1 beta1 beta1 alpha1; ...
                  alpha alpha beta gamma; ...
                  alpha alpha gamma beta; ...
                  alpha beta alpha gamma; ...
                  alpha beta gamma alpha; ...
                  alpha gamma beta alpha; ...
                  alpha gamma alpha beta; ...
                  beta alpha alpha gamma; ...
                  beta alpha gamma alpha; ...
                  beta gamma alpha alpha; ...
                  gamma alpha alpha beta; ...
                  gamma alpha beta alpha; ...
                  gamma beta alpha alpha];                  
        weight = [0.1884185567365411, ...
                  w1, w1, w1, w1, ...
                  w2, w2, w2, w2, w2, w2, ...
                  w2, w2, w2, w2, w2, w2];
end
%% Verification
% The order of the quadrature rule is verified by the function
% verifyquadpts. See <matlab:ifem('verifyquadpts') verifyquadpts>.
